%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% This program replicates Figure 1 in the paper "Optimal 
%%% Portfolio Choice with Fat Tails and Parameter Uncertainty",
%%% Journal of Financial and Quantitative Analysis (2024),
%%% Raymond Kan & Nathan Lassance 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
close all

%%% Exponential integral
Expfun = @(n,x) integral(@(t)exp(-x.*t)./(t.^n),1,inf);

%%% Initialize vectors
rhovec = 0:0.01:1;
etanu4 = zeros(length(rhovec),1);
etanu8 = zeros(length(rhovec),1);
varphinu4 = zeros(length(rhovec),1);
varphinu8 = zeros(length(rhovec),1);

%%% for loop
j = 1;
for rho = rhovec
    %%% nu = 4
    nu = 4;
    if rho == 0
        etanu4(j,1) = 1;
        varphinu4(j,1) = 1;
    elseif rho == 1
        etanu4(j,1) = nu/(nu-2);
        varphinu4(j,1) = nu^2/(nu-2)^2;
    else
        y = @(x) (nu-2).*rho.*x./(2*(1-rho));
        fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
        etanu4(j,1) = fzero(fun,0.5*(1+nu/(nu-2)));
        varphinu4(j,1) = 2*etanu4(j,1)^2*(1-rho)/(nu-etanu4(j,1)*(nu-2));
    end
    %%% nu = 8
    nu = 8;
    if rho == 0
        etanu8(j,1) = 1;
        varphinu8(j,1) = 1;
    elseif rho == 1
        etanu8(j,1) = nu/(nu-2);
        varphinu8(j,1) = nu^2/(nu-2)^2;
    else
        y = @(x) (nu-2).*rho.*x./(2*(1-rho));
        fun = @(x) y(x).*exp(y(x)).*Expfun(nu/2,y(x))-rho;
        etanu8(j,1) = fzero(fun,0.5*(1+nu/(nu-2)));
        varphinu8(j,1) = 2*etanu8(j,1)^2*(1-rho)/(nu-etanu8(j,1)*(nu-2));
    end
    j = j+1;
end

%%% Create Figure 1
fig = figure;
plot(rhovec,etanu4,'b','LineWidth',2);
hold on
plot(rhovec,etanu8,'-.b','LineWidth',2);
plot(rhovec,varphinu4,'--r','LineWidth',2);
plot(rhovec,varphinu8,':r','LineWidth',2);
hold off
legend('$\eta$ $(\nu=4)$','$\eta$ $(\nu=8)$','$\varphi$ $(\nu=4)$',...
       '$\varphi$ $(\nu=8)$','Location','northwest','interpreter','latex');
set(gca,'TickLabelInterpreter','latex');
grid on
xlabel('$\rho$','interpreter','latex');
xlim([0 1]);
ylim([1 4]);
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([1 1.5 2 2.5 3 3.5 4]);
fontsize(fig,20,"points");
